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Abstract 

We consider the numerical discretization of the time-domain Maxwell's 
equations with an energy-conserving discontinuous Galerkin finite element 
formulation. This particular formulation allows for higher order approxi- 
mations of the electric and magnetic field. Special emphasis is placed on 
an efficient implementation which is achieved by taking advantage of re- 
currence properties and the tensor-product structure of the chosen shape 
functions. These recurrences have been derived symbolically with com- 
puter algebra methods reminiscent of the holonomic systems approach. 



1 Introduction 

This paper is dedicated to a successful cooperation between symbolic compu- 
tation and numerical analysis. The goal is to simulate the propagation of elec- 
tromagnetic waves using finite element methods (FEM). Such simulations play 
an important role for constructing antennas, electric circuit boards, bodyworks, 
and many other devices where electromagnetic radiation is involved. The nu- 
merical simulation of such physical phenomena helps to optimize the shape of 
components and saves the engineer from doing a long and expensive series of 
experiments. 

*This article is part of the volume U. Langer and P. Paule (eds.) Numerical and Sym- 
bolic Scientific Computing: Progress and Prospects in the series Texts & Monographs in 
Symbolic Computation, ISBN 978-3-7091-0793-5. The original publication is available at 
www.springerlink.com, DOI 10.1007/978-3-7091-0794-2.6. 

"supported by the Austrian Science Fund (FWF): SFB F013 and P20162-N18, and partially 
by NFS-DMS 0070567 as a postdoctoral fellow. 
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Finite clement methods serve to approximate the solution of partial differen- 
tial equations on a given domain C R d subject to certain constraints (e.g., 
boundary conditions). The domain f2 is partitioned into small elements (typi- 
cally triangles or tetrahedra) and the solution is approximated on each element 
by means of certain shape functions. In our application we deal with Maxwell's 
equations which relate the magnetic and the electric field. In Section 2 we de- 
scribe how the problem can be discretized using FEM and in Section 3 we give 
the details concerning an efficient implementation. 

An important ingredient for the fast execution of some operations in the FEM 
are certain difference-differential relations that were derived with computer alge- 
bra methods. The methods that we employ, originate in Zcilberger's holonomic 
systems approach [13, 3, 10] whose basic idea is to define functions and se- 
quences in terms of differential equations and recurrence equations plus initial 
values (these equations have to be linear with polynomial coefficients). Luckily 
the shape functions used in the chosen FEM discretization fit into the holonomic 
framework since they are defined in terms of orthogonal polynomials. Section 4 
explains how the desired relations have been computed. 



2 FEM formulation of Maxwell's equations 

In order to describe electromagnetic wave propagation problems, we consider 
the loss-free time-domain Maxwell's equations 

e— = curl if, 
at 

li— = -cw\E, 

subject to appropriate initial and boundary conditions. Here E = E(x,t) 
denotes the electric and H = H(x, t) the magnetic field strength (with x = 
(xi,X2,X3) the space variables and t the time), and e and fi > arc the per- 
mittivity and the permeability, respectively. When discretizing these equations 
with the finite element method, we go over to a weak formulation by multiplying 
both equations with test functions e(x) and h(x) and integrating over the whole 
domain il C R 3 . The solution of the Maxwell's equations then has to fulfill the 
conditions 

d 

— (eE,e) n = (cm\H,e) Q , 

8 (1) 

— (pH,h) n = -(cmlE,h)n 

for all test functions e and h, where (•, -)q is the short notation for the L 2 (f2) 
inner product (a, 6)q — J Q abdx. Then we replace both the magnetic and 
electric field as well as the test functions by finite-dimensional approximations 
on a triangulation 7~h of the domain f2. Herein h denotes some characteristic 
length of the elements in Th (not to be confused with the test function h). 

Conforming finite elements ensure that the finite-dimensional approximations 
are within a space which is appropriate for the partial differential equations 
under consideration. For Maxwell's equations this space is if (curl, fi) which de- 
mands tangential components to be continuous across element interfaces. The 
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discontinuous Galcrkin finite element method (DG) neglects this conformity 
condition when building up a discrete basis for the approximation, but instead 
has to incorporate stabilization terms to achieve a consistent and stable formu- 
lation. This is normally done by applying integration by parts and replacing 
fluxes at element boundaries with numerical fluxes [1, 11, 8, 7]. The latter ap- 
proach has the major advantage that the mass matrices M e and M^, i.e., the 
matrices that arise when discretizing (sE, e)n and ([iH y h)ci, respectively, are 
block-diagonal which makes the application of their inverses computationally 
more efficient. 

We consider the approximation space 

V h k = {ve (L 2 (n)) 3 : v\ T G (V k (T)) 3 VT e %} 

that consists of functions which are piecewise polynomial up to degree k. By 
integration by parts of (1) on each element T € Th, and by adding a consistent 
stabilization term on all element boundaries we get (again for all test functions e 
and h) 

^ ]T (eE, e) T = (( H > curl e )t + ( H * x e) 9T ) , 

TeT h Ten 

£ ]T ^H, h) T = (- ( curl E > + (E*-E,hx u) dT ) , 

Ten Ten 

where v denotes the outer normal on each element boundary and H* , E* are the 
numerical fluxes. The properties of different DG formulations mainly depend 
on the choice of the numerical fluxes. As all derivatives are now shifted to the 
electric field E and the according test functions e, it is reasonable to approximate 
the electric field of one degree higher than the magnetic field. So we choose the 
approximation spaces V^ +1 for E and e and V* for H and h. 



2.1 Numerical flux 

Several choices for the numerical flux are used in practice. Our goal here is 
to derive a numerical flux which ensures that the numerical approximation ful- 
fills the following two important properties which are already fulfilled on the 
continuous level: 

1. conservation of the energy ^(eE,E)q + ^(^,H,H)q 

2. non-existence of spurious modes 

On the one hand using dissipative fluxes avoids spurious modes and is often 
used, but as it introduces dissipation, the energy of the system is not conserved. 
On the other hand the standard approach for energy conserving methods is the 
so called central flux. Its mayor disadvantage is, that it introduces non-physical 
modes, spurious modes. 

Nevertheless we start with this approach to derive the stabilized central flux 
formulation which gets rid of both problems. A more extensive discussion of 
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numerical fluxes (including the stabilized central flux) for Maxwell's equations 
can be found in [6, 8.2]. 

The central flux takes the averaged values of neighboring elements for the nu- 
merical flux, i.e., H* = fHj and E* = fEj with -{{•]} denoting the averaging 
operator, and ends up with a semi-discrete system of the form 



d (M F \ (E\ ( -C%\ (E 



(2) 



where Ch denotes the discrete curl operator stemming from the central flux 
formulation. The matrix on the left side is symmetric and positive definite 
whereas the matrix on the right side is antisymmetric. Then the evolution 
matrix for the modified unknowns (Mi E 1 H) T is also antisymmetric and 
thus the proposed energy is conserved. Nevertheless this matrix has a lot of 
eigenvalues close to zero which correspond to the discretization, but not to 
the physical behavior of the system. To motivate the modification which will 
stabilize the formulation, let us have a brief look at the problem in frequency 
domain, i.e., for time-harmonic electric and magnetic fields. Then the discrete 
problem in frequency domain reads (with frequency lo): 

= (ilo) 2 (M £ E, e) + (M^ChE, C h e). (3) 

The problem with non-physical zero eigenvalues now manifests in (ChE, Che) 
being only positive semidefinite. We overcome this issue by adding a stabiliza- 
tion bilinear form S(E,e) to (3) as proposed in [6]. 



S{E,e) := ^(M^.W^V 



a 

Fer h 



with a > 0, where Th is the union of all element boundaries and [•] denotes 
the jump operator, i.e., the difference between values of adjacent elements. This 
stabilization bilinearform eliminates the nontrivial kernel of Ch and is consistent 
as \E\ x v is zero for the exact solution. Before we can translate the formulation 
back to the time domain, we introduce a new variable which is defined as 

h f , = (jgl x v ) a 

The new unknown H F is also piecewise polynomial on each face. 

If we go back to the time-domain formulation we end up with the following 
formulation (note that relations between [•] and -{{•]}• were used): 

^ ]T (sE, e) T - ((#, curl e) T + ({H} x u, e) dT ) + 

TeT h TeT h 

Fer h 



^ Y (l*H, h) T = Y ( - curl E )t + {\{E\xv, h) dT ) , 
TeT h TeTh 



dt J-L h 

FeT h Fer h 
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For p-robust behavior a should scale with p 2 , where p is the polynomial degree. 
This is motivated by the symmetric interior penalty method for elliptic equations 
(see e.g. [1]) where a scaling of a with p 2 in the bilinearform S is necessary for 
stability to dominate over some terms stemming from inverse inequalities which 
scale with p 2 (see also [7]). 

We again achieve a system of the form (2) where the vector H now consists 
of element and face unknowns and the matrix representing the discrete curl 
operator is the stabilized central flux curl operator now. Thus we conclude 
that the method now conserves energy, and spurious modes, introduced by the 
central flux, vanish. 



2.2 Numerical Examples (Spherical Vacuum Resonator) 

We consider a spherical domain £1 := {x £ R 3 : ||x||2 < 1} and the frequency do- 
main formulation of the Maxwell's equations subject to perfect electrical bound- 
ary conditions 



iujeE — curl H, 
iujfiH = — curl_E, 
Exv = 



on 
on 



n. 



To demonstrate the opportunities of higher order discretizations we consider 
a coarse mesh consisting of 30 elements and increase the polynomial degree 
to increase the spatial resolution. We are interested in the error of the eight 
smallest resonance frequencies. Therefore we compare the eigenvalues of the 
numerical discretization with those of a reference solution. In Figure 1 we 
observe the expected exponential convergence of the method. 



Testcase 2 - relative error of resonant frequencies 
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Figure 1: Convergence of the resonance frequencies after p-refinement 
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3 Computational aspects 



As the spatial discretization conserves energy, we consider symplectic time in- 
tegration methods which conserve the energy on a time-discrete level. The 
simplest one is the symplectic Euler method which discretizes the semi-discrete 
system (2) in the following way: 

H n+l = jjn + At M -l ChE n 
E n+1 = E n _ At M - X ClH n+1 

with the stability condition 

te<2( K p{M- 1 *C h M- 1 CZMP)} 

The matrix M M 2 C^M^C^M^ 2 is symmetric and the spectral radius p can 
be estimated once by an iterative method like the power iteration. 
When shifting the electric or the magnetic field by a half time-step we can re- 
construct the well-known leap frog method. Nevertheless for our considerations 
it is less important which time integration scheme is used as long as it is explicit. 
The matrix multiplications with Ch and Cj^ (see Section 3.2) as well as with 
M~ x and M~ x (see Section 3.3) decide about the computational efficiency of 
an implementation. 

The advantage of discontinuous Galerkin methods becomes evident now. The 
mass matrices can be inverted in an element by element fashion and also the 
discrete curl operations only need information of (element-)local and adjacent 
degrees of freedom, which allows for straightforward parallelization. Element 
matrices such as mass matrices and the discrete curl operation can be stored 
once and applied at each time step. This is how far one comes just because of 
the formulation itself. 

With appropriate choices for the local shape functions we can use advanced tech- 
niques to execute those operations with a lower complexity than local matrix- 
vector multiplications. Furthermore we don't even have to store the element 
matrices, s.t. the techniques presented below are also much more memory- 
efficient. 

The following ingredients are essential for the techniques proposed below, which 
enhance the implementation of the DG method: 

1. Definition of an L 2 -orthogonal basis of polynomial shape functions in 
tensor-product form 1 on a reference clement T 

2. Use of curl-conforming (covariant) transformation for evaluations on the 
physical element T 

3. Use of recurrences for the polynomial shape functions to evaluate gradients 
and curls 

4. Use of tensor-product structure to evaluate traces 2 

1 these arc polynomials which arc products of univariate polynomials 

2 values at a boundary 
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3.1 Local shape functions 



For stability and fast computability we choose the ^-orthogonal Dubiner basis 
[5, 9]. Here, the basis functions on the reference element are constructed in a 
tensor-product form of Jacobi polynomials P^ a '^ for each spatial component 
(note that the Legendre polynomials P, = P> ' ' are just a special case). For 
example, on the reference triangle spanned by the points (0, 0), (1, 0) and (0, 1) 
the shape functions take the form 

<p id (x,v) = Pi - l) • (1 - xf ■ Pf l+lfi \2x - 1). (4) 

They are orthogonal on the reference triangle, and gradients can be evaluated 
by means of recurrence relations as demonstrated in Section 3.2.2. Due to the 
tensor-product form traces can be evaluated very fast, see Section 3.2.3. 



3.2 Discrete curl operations 

At each time step we have to evaluate terms like (H, curle)r on each element T 
and {^H^ x u, [e])i? on each face F. Similar expressions have to be evaluated 
for the electric field E. 



3.2.1 Covariant transformation 



Let $ : f — > T be a diffcomorphic mapping from the reference element to some 
physical element T. Then the covariant transformation of a function u defined 
on the reference element T is 



u := (-F -1 ) T uo$- 



with 



F = V$. 



If we define the shape functions on the mapped elements as the covariant trans- 
formed shape functions on the reference element, then the tangential compo- 
nent on the mapped element depends only on the tangential component of 
the reference element. The transformation is called curl-conforming as it en- 
sures that for any function u G i?(curl, f2) the covariant transformed func- 
tion u lies in if (curl, f2). Furthermore it preserves certain integrals, s.t. the 
following relations hold for the covariant transformations H , e G H (curl, T) of 
H,e e ff(curl,T): 



L 

/< 

JdT 



H curl e da; 



curl e Ax 



(H x v)eds 



/ (H x v)eds 

Jdt 



This means that the integrals of these forms appearing in the formulation are 
independent of the geometry of the particular elements. The matrices can be 
computed once on the reference clement. This trick was published in [4]. 



7 



3.2.2 Evaluating gradients 



For computing curls it is sufficient to evaluate gradients, since the curl is a 
certain linear combination of derivatives. We write the corresponding function E 
in modal representation, i.e., 

E = ^2a a Lp ai a a e R 3 , 

where the sum ranges over the finite collection of (scalar) shape functions defined 
on the reference element (in 2D the multi-index a is and in 3D a = (i,j,k)). 
With the use of the covariant transformation, we just have to consider the 
integral on the reference element T: 

hcmlEdx. 

The idea is now to take advantage of recurrence relations between derivatives 
of Jacobi polynomials and Jacobi polynomials itself. We aim for an operation 
which gives the coefficients b a € R 3 representing the gradient 

VE = ^2b a tp a . 

a 

Then L 2 -orthogonality can be used to evaluate the complete integral very fast. 

For ease of presentation let's consider the far more easy case of evaluating the 
derivative of a scalar one-dimensional function v(x) — Y^ii=o v i p i( x )i Vi £ R 
given in a modal basis of Legendre polynomials Pi, which fulfill the relation 

P' i+1 {x)=P' i _ 1 {x) + {2i + l)P i {x). (5) 

Then the problem is to find the modal representation of 

n n—1 

v\x) = Y,ViPl{x) = Y, w i p i( x )- 

i=0 j=0 

Let's show the first step, i.e., how we get the highest order coefficient w n -\\ 

n n — 1 

i=0 i=0 
n-1 

= V ^ X ) + V nK-2( X ) + Vn{2n - 1)P„_! 

i=0 
n-1 

= J! ViPi{x) + w„-iP„-i(x) 

i=0 

where we used the recurrence relation (5) for -P^(x) and thus get w„-i = v n (2n— 
1). For the remaining polynomial ^2™=q ViPl(x) of degree n— 1 we can apply the 
same procedure to get w„-2- This can be continued until also wo and thereby 
the complete polynomial representation ^2™=q WiPi(x) of v'{x) is determined. 
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An efficient C++ implementation of this procedure was achieved by template 
meta-programming, where the compiler can generate optimized code for all ele- 
ments up to an a priori chosen maximal polynomial order. 

The same basically also works in three dimensions with Jacobi polynomials, but 
the relations are far more complicated, see Section 4, and need 3 nested loops. 

The overall costs for the evaluation of the element curl integral scales linearly 
with the number of unknowns N on one element which is much better than the 
matrix- vector multiplication which already has complexity 0(N 2 ). 

3.2.3 Evaluating traces 

The boundary integrals that have to be evaluated can make use of the tensor- 
product form to evaluate traces. Again we don't want those traces to be evalu- 
ated pointwise but in a modal sense and recurrences for the Jacobi polynomials 
make the transformation from volume element shape functions to face shape 
functions with O(N) operations possible. The procedure therefore is similar to 
the evaluation of the gradient in the previous section. 

3.3 Mass matrix operations 

So far we dealt only with the discrete curl operations. So the only thing that 
is left to talk about is the application of the inverse mass matrices. Due to the 
covariant transformation we have 

{{M e ) a .f)) Lm = J e((p a e[)(ip /3 e m )dx 

= J^\& C t{F)\e{^ a ef)F- 1 {F' 1 ) T ^pe m )dx (6) 

with ip a denoting the scalar-valued shape functions and e„ the n-th unit vector. 
Note also the block structure of M £ that is indicated by the above notation. 
In some FEM applications, symbolic methods related to those described in 
Section 4, can be used to prove the sparseness of the corresponding system 
matrix, see [12]. 



3.3.1 Flat elements 

Let's assume the material parameters e and \x are piecewise constant and the 
elements are flat, i.e., V$ = F = const on each element. Then the integral (6) 
simplifies to 

/ e((p a e[) ((p e m )dx = | det(f )| e (F" 1 (F" 1 ) T ) Zjm / tp a <ppdx 
Jt Jf 

and as Jf <~p a '~pp da; = 8 a ^ the matrix is (3 x 3)-block-diagonal and the inversion 
is trivial. The computational effort is obviously of order O(N) where N is the 
number of unknowns. 
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3.3.2 Curved elements 



If we consider curved elements or non-constant material parameters e and /U, 
the approach has to be modified as the mass matrix arising from (6) may be 
fully occupied. Let's go a step back and consider a similar scalar problem 3 with 
a non-constant coefficient e: 



Given: f( v )~ J fvdx 

Find: u, s.t. / suvdx = / fvdx 
Jt Jt 

We now transform back to the reference element T and get 

/ suvdx — / \ det(F)\ euv dx = / uvdx 
Jt Jt Jt 

where v = \ det(F) \ ev. If we now approximate v with the same basis we used 
for v before, the mass matrix is diagonal again. Nevertheless the evaluation of 
the functional f(v) has to be transformed as well: 

fvdx — f — — \- —— fvdx = f -fv dx 
J T \det(F)\e J J f e J 

To evaluate the last term we will use numerical integration. But as (in our 
application) / is not given pointwise, but in a modal sense, we have to calculate 
a pointwise representation for the numerical integration of J T fv dx first: 

Given: f(v) = [ fv dx = [ | dct(F)| fv dx 
Jt Jt 

Find: f u s.t. / fvdx = ^ | det(F)\(xi)fjUJiv(xi) 

Then we can divide (on each integration point) by e and with those new coeffi- 
cients we can, by numerical integration, get a good approximation to f T ^fvdx. 
The "reverse numerical integration" and the numerical integration used here can 
be accelerated by the use of the sum factorization technique. Doing so the com- 
plexity of both "reverse numerical integration" and the numerical integration 
is 0(p ), where p is the polynomial degree. Note that the approximate inverse 
M" 1 obtained by this method is still symmetric and positive definite. 



3.4 Overall computational effort 

In the previous sections we saw that the overall computational effort scales 
linearly with the degrees of freedom N as long as the elements are flat and 
coefficients are piecewise constant. Even for curved elements (and variable co- 
efficients) the computational effort is only of order O(N^). Furthermore no 
element matrices have to be stored. Only the geometric transformations and 
the local topology have to be kept in the memory. 

Extensions to 3D arc straightforward 
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3.5 Timings 



Let's also state some exemplary numbers that were achieved for this method 
and its implementation on an Intel Xeon CPU 5160 at 3.00 GHz (64 bit) (single 
core) for a tetrahedral mesh with 2078 elements. The costs for one step of the 
symplectic Euler method per 6 scalar degrees of freedom are listed in Table 1. 



order p 


time [usee] 


order p 


time [usee] 


1 


0.61 


1 


4.89 


2 


0.58 


2 


2.54 


3 


0.71 


3 


1.93 


4 


0.79 


4 


1.79 


5 


1.16 


5 


2.06 


6 


1.24 


6 


2.17 


7 


1.32 


7 


2.33 


8 


1.53 


8 


2.67 


9 


1.66 


9 


2.88 


10 


1.74 


10 


3.04 



Table 1: Timings for flat elements (left), using 0(1) floating point operations 
per dof and curved elements (right) using 0(p)) floating point operations per 
dof. 



4 Symbolic derivation of relations 

In this section we want to describe the symbolic methods that were employed 
for finding the desired relations for the polynomial shape functions. These 
relations allow for efficient computation of the discrete curl operations and 
traces as described in Section 3.2. They have been computed by following 
the holonomic systems approach [13, 3, 10], which works for all functions that 
satisfy sufficiently many linear differential equations or recurrences or mixed 
ones; these relations have to have polynomial coefficients. A large class of func- 
tions (like rational or algebraic functions, exponentials, logarithms, and some 
of the trigonometric functions) as well as a multitude of special functions is 
covered by this framework. Part of it are algorithms for the "basic arithmetic" 
(that we will refer to as "closure properties"), i.e., given two implicit descrip- 
tions for functions / and g, respectively, we can compute such descriptions 
for / + <?, fg, and for functions obtained by certain substitutions into / or 
g. All computations in this section have been performed in Mathematica us- 
ing our package HolonomicFunctions (it is freely available from the website 
http : / / www. rise . uni-linz . ac. at / research/combinat /software / ) . 

4.1 Introductory example 

For demonstration purposes we show how to derive automatically the rewriting 
formula (5) for Legendre polynomials P n (x) . It is well known that these orthog- 
onal polynomials satisfy some linear relations, e.g., the second order differential 
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equation 

(x 2 - l)P^(x) + 2xP' n {x) - n(n + l)P n (x) = 
or the three term recurrence 

(n + 2)P n+2 {x) - (2n + 3)xP n+1 (x) + (n+ l)P n {x) = 0. 

We will represent such linear relations in the convenient operator notation, using 
the symbols D x for the partial derivative with respect to x, and 5^ for denoting 
the shift operator with respect to n. Then the two relations above are written 
as 

(x 2 - l)D 2 + 2xD x + (-n 2 - n) 

and 

(n + 2)S 2 + {-2nx - 3x)Sn + (n + 1), 

respectively, and we identify operators and relations with each other. The op- 
erators can be regarded as elements of a (noncommutative) polynomial ring in 
Sn and D x with coefficients being rational functions in Q(n, x). We can obtain 
additional relations for P n (x) by combining the given relations linearly, or by 
shifting and differentiating them. In the operator setting these operations cor- 
respond to addition and multiplication (from the left) and we can refer to the 
set of all operators obtained in this way as the annihilating left ideal generated 
by the initially given operators. In the following we will represent annihilating 
ideals by means of their Grobner bases; these are special sets of generators that 
allow for deciding the ideal membership problem (i.e., the question whether 
some relation is indeed valid for the function under consideration) and for ob- 
taining unique representatives of the residue classes modulo the ideal (see [2]). 
All algorithms mentioned below will require Grobner bases as input. A Grobner 
basis of the annihilating ideal of the Legendre polynomials is given by 

G={(n+ l)Sn + (1 - x 2 )D x + {-nx - x), {x 2 - l)D 2 + 2xD x + (-n 2 - n)}. 

Our main task will be to find elements with certain properties in an annihilating 
ideal; this can be done via an ansatz as we demonstrate now. The relation (5) 
that we are going to recover connects P' n+2 {x), Pn(x), and P n+ \(x), and its 
coefficients are free of x. These facts translate to an ansatz operator of the form 

A = Cl (n)D x S 2 + c 2 (n)D x + c 3 (n)^ 

where the coefficients c, are rational functions in Q(n), and hence free of x as 
required. We have to determine the Cj such that the operator A is an clement 
of the left ideal / generated by G, so that A(P n (x)) = 0. For this purpose we 
use the Grobner basis G to compute the unique representation of the residue 
class of A modulo I (it is achieved by reduction) . We have A € I if and only if 
the residue class is represented by the zero operator and hence we can equate 
all its coefficients to zero, obtaining the following two equations 

Cl {2nx 2 + 3x 2 - n- 2) + c 2 (n + 1) + c 3 (x 2 - 1) = 0, 
c 1 (n+l)(2n + 3)x + c 3 (n + l)x = 0. 
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Note that in these equations the variable x occurs, since it is contained in the 
coefficients of G. We get a solution that is free of x by performing a coefficient 
comparison with respect to this variable. This yields in the end the linear system 

-n - 2 n+l 

2n + 3 J I I c 2 I = 




^(n+l)(2n + 3) 
whose solution is 

Ci = — 1, c 2 = 1, c 3 = 2n + 3, 
and this gives rise to the desired relation. 

Now what do we do if we don't know the exact shape of the ansatz as given 
here by A? Then we have to include all possible monomials D*S^ up to some 
total degree into our ansatz. Looping over the degree, we will finally find the 
relation, but the effort can be tremendous. Therefore, as a preprocessing step, 
we determine the shape of the ansatz by modular computations. This means 
plugging in concrete values for some of the variables and reducing all integers 
in the coefficients modulo some prime. These techniques have been described 
in detail in [10] and they are crucial for getting results in a reasonable time. 

All these steps have been implemented in the package HolonomicFunctions and 
it computes the relation (5) immediately: 
in[i] := « HolonomicFunctions. m 

HolonomicFunctions package by Christoph Koutschan, RISC-Linz 
Version 1.3 (25.01.2010) 

Type ?HolonomicFunctions for help 



inp]:= FindRelation[Annihilator[LegendreP[n, x]], Eliminate — >■ x\ 

ou,[2]= {S*D X + (-2n - 3)& - D x } 



4.2 Relations for the shape functions 

A core functionality of our package HolonomicFunctions [10] is to execute 
closure property algorithms (e.g., for addition, multiplication, and substitution) 
on functions represented by their annihilating ideals. We can now use these 
algorithms to obtain annihilating ideals for the shape functions ip, since their 
definition in terms of Jacobi and Legendre polynomials involves just the above 
mentioned operations. 



4.2.1 The 2D case 

We first consider triangular finite elements in two dimensions. For these, the 
shape functions are defined as in (4). Analogously to the one-dimensional ex- 
ample in Section 3.2.2 we want to express the partial derivatives (with respect 
to x and y, respectively) in terms of the original shape functions. So the goal is 
to find relations (free of x and y) that connect the partial derivatives with the 
original function. More concretely, we are looking for a relation that allows to 
express some linear combination of shifts of ^ipij(x, y) as a linear combination 
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of shifts of ipi.j(x,y) (and similarly for y). This corresponds to an operator of 
the form 

E c hm!n (i,j)D x srs J n + c o,mAhj)srs? (7) 

(m,n)GN 2 (m,n)eM ! 

where the yet unknown coefficients Cd, m ,n € Q(z, j) do not depend on x and y, 
and the sums have finite support. 

Since we have to find such a relation in the annihilating ideal for <pij(x,y), 
it is natural to start by computing a Grdbner basis for this ideal. The pack- 
age HolonomicFunctions provides a command Annihilator that analyzes a 
given mathematical expression and performs the necessary closure properties 
for obtaining its annihilating ideal. So in our example we can just type 

inp] = ann = Annihilator[(l — x)"i * LegendreP[i, 2y/(l — as) — 1] * 

JacobiP[j,2i + 1,0,2a;- l],{S[i], S[j], Der[a;], Der[y]}]; 

and after a second we have the result (which is already respectable in size, 
namely 340kB, corresponding to about 10 pages of output). 

Having implemented noncommutative Grobner bases, our first attempt was to 
use them for eliminating the variables x and y. But it soon turned out that 
this attempt did not produce optimal results, and in addition the computations 
were very time-consuming. Therefore we came up with the ansatz described in 
Section 4.1. We use it now to compute the desired relations (both computations 
take less than a minute): 

in[4] = FindRelation[ann, Eliminate — >• {x, y}, Pattern — >• {_, _, | 1, 0}] 
// Factor 

out[4]= {(2i +j + 5)(2i + 2j + 5)SiS 3 2 D x + (j + 3)(2i + 2j + 5)S J 3 D X + 
2{2i + 3){i + j + 3)S t S ] D x - 2(2i + + j + 3)S 2 D X - 
2{i + j + 3)(2i + 2j + 5)(2i + 2j + 7)SiSj - (j + l)(2i + 2j + 7)SiD x - 
2{i +j + 3)(2i + 2j + 5)(2i + 2j + 7)S- - (2i +j + 3)(2i + 2j + 7)S j D x ) 

in[5]:= FindRelation [ann, Eliminate — y {cc, y}, Pattern — >• {_, _, 0,0 | 1}] 
// Factor 

ou,[ 5 ]= {(2i + j + 6)(2i+j + 7){2i + 2j + 7)S l 2 S J 2 D y - (j + 3)(j + 4)(2t + 2j + 7)S?D y - 
4(j + 2)(i + j + 4)(2i + j + 6)S 2 S 3 D y + 4(j + 3)(t + j + 4)(2i + j + 5)S J 3 D y + 
{j + 1) (j + 2) (2t + 2j + 9)S 2 D y - 4(2t + 3) (t + j + 4) (2t + 2j + 7) (2i + 2j + 9)SiS 2 - 
(2i + j + 4) (2i +j + 5) (2i + 2j + Q)S?Dy} 

Here the option Pattern specifies the admissible exponents for the operators, 
e.g., in the first case we allow any exponent for the shift operators, whereas D x 
may occur with power at most 1 only, and D y must not appear at all in the 
result. 

4.2.2 The 3D case 

When dealing with tetrahedra in three dimensions, the shape functions are 
denoted by <Pij,k{x, y, z) and are defined by 

(1 - x y)'(l - xyP^^ l) /f <+1 '°> (£_ _ !) P f +2 J+2 ,o )(2x _ 1} _ 
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Again they have the nice property of being L 2 -orthogonal on the reference tetra- 
hedron 

T = {(x, y,z)eB 3 \x>0Ay>0Az>0Ax + y + z<l}. 

Computing an annihilating ideal for <Pij,k(%, V, z) is already much more involved 
than in the 2D case: 

in[e] := phi = (1 — x — y)"i (1 — x)"j LegendreP[i, 2z/(l — x — y) — 1] 

JacobiP|>', 2i + 1, 0, 2y/(l - x) - 1] JacobiP[fe, 2i + 2j + 2, 0, 2x - 1]; 
in[7]= Timing[ann = Annihilator[phi, {Der[cc], S[i], S[j], S[fe]}];] 

ou,[7]= {359.686, Null} 

The Grobner basis for this annihilating ideal is about 117MB in size (correspond- 
ing to several thousand of printed pages). Note also that it is more efficient to 
consider only one derivation operator, and compute annihilating ideals for each 
of the cases and ^ separately (this applies to the 2D case, too). 

In principle, the desired relations for the 3D case can be found in the same way 
as for two dimensions. As described in Section 4.1 we find by means of modular 
computations that the ansatz (for the case ^) contains the 16 monomials 

SiSjSk-D x , SiS^D xl S^S^D X1 SjSj?D x , SiSjSkD x , SiS^D xl S^SkD x , SjS^D x , 
SiSjSk, SiSjD x , StS k , SiSkD x , Sj Sk, S) D x , SjS k , SjSkD x . 

However, in order to compute the corresponding coefficients, we did not succeed 
with the standard approach used in Section 4.2.1. Instead, we had to employ 
modular techniques again for many interpolation points, and then interpolate 
and reconstruct the solution. 

5 Conclusion 

We have presented an efficient implementation for solving the time-domain 
Maxwell's equations with a finite element method that uses discontinuous Ga- 
lerkin elements. Besides many other optimizations that speed up the whole 
simulation, the usage of certain recurrence relations for the shape functions al- 
lows for a fast evaluation of gradients and traces. These relations have been 
derived symbolically with computer algebra methods. 

It is widely believed that the mathematical subjects "numerical analysis" and 
"symbolic computation" do not have much in common, or even that they are 
kind of orthogonal. Experts from both areas can barely communicate with 
each other unless they don't talk about work. It was the great merit of the 
project SFB F013 "Numerical and Symbolic Scientific Computing" that had 
been established in 1998 at the Johannes Kepler University of Linz, Austria, to 
bring together these two communities to identify potential collaborations. We 
consider our results as a perfect example for such a fruitful cooperation. 
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